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I ■  INTRODUCTION 


The  damaging  effect  of  low  velocity  impact  upon  laminated  composite 
structures  has  been  recognized  by  the  engineering  community  as  an  area  of 
primary  concern.  This  type  of  loading  can  easily  be  envisioned  as  a 
workmans  tool  being  dropped  on  a  structure  during  assembly,  or  runway  debris 
thrown  against  a  composite  airframe  during  takeoff  or  landing.  The  damage 
produced  by  these  loads  may  consist  of  matrix  cracking,  delaminations  and 
even  fiber  failures,  depending  on  the  mass  and  velocity  of  the  impactor. 
Often,  this  damage  is  invisible  to  the  naked  eye,  and  may  go  unnoticed  even 
though  the  strength  of  the  structure  has  been  seriously  degraded. 

The  primary  objective  of  Phase  1  of  this  program  was  the  development  of 
an  analytical  model  for  studying  the  transient  dynamic  response  of  laminated 
composites  during  this  type  of  impact  event.  This  model  utilized  a  plate 
bending  finite  element  to  characterize  the  laminate  and  a  lumped  mass  and 
nonlinear  contact  force  spring  to  characterize  the  impactor.  The 
discretized  equations  of  motion  were  generated,  and  integrated  numerically 
to  yield  the  displacement  history  of  the  system.  The  displacement  history 
could  then  be  used  to  determine  the  strain,  and  ultimately  stress  histories 
for  the  laminate.  These  stresses  and/or  strains  were  then  used  to  predict 
the  onset  of  damage  in  the  plate  (Ref.  1). 

The  Phase  I  model  predicts  stresses  and  strains  which  are  consistent 
with  the  plate  bending  assumptions.  That  is,  extensional  strains  vary 
linearly  through- the -thickness  of  the  plate,  and  shear  strains  vary 
parabolically .  The  distributions  predicted  by  the  model  are  then  a  good 
representation  of  the  far  field  solution,  but  do  not  adequately  describe  the 
complex  fields  that  may  exist  in  close  proximity  to  the  impactor. 

Accurate  characterization  of  these  stress  and  strain  fields  is 
important  since  they  may  result  in  peak  values  and  distributions  that  are 
substantially  different  from  those  predicted  by  plate  bending  theory. 
Development  of  a  more  refined  model  for  predicting  the  three-dimensional 
state  of  stress  in  the  vicinity  of  the  impactor  has  been  the  focus  of  the 
Phase  II  work. 

This  work  has  concentrated  on  three  topics;  development  of  a  three- 
dimensional  laminated  solid  finite  element,  optimization  of  the  algorithms 
used  to  generate  and  solve  the  equations  of  motion,  and  validation  of  the 
resulting  analytical  model. 
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II ■  ANALYTICAL  HODEL 


1.  OVERVIEW 

The  Intent  of  the  analytical  model  is  to  provide  an  accurate  tool  for 
predicting  the  initiation  of  damage  in  a  laminated  composite  due  to  low- 
velocity  impact.  The  definition  of  low-velocity  impact  is  debatable,  but 
for  this  program  it  will  be  interpreted  as  an  impact  event  which  produces 
dynamic  effects  that  are  dominated  by  flexural  waves. 

To  allow  generality  in  defining  parameters  such  as  plate  geometry  and 
boundary  conditions,  finite  elements  are  used  to  model  the  laminate.  The 
impactor  is  assumed  to  be  very  stiff  relative  to  the  plate  allowing  it  to  be 
treated  as  a  lumped  mass  coupled  to  the  plate  through  some  suitable  contact 
lav . 

It  was  assumed  that  all  material  non-linearities  would  only  be  apparent 
in  a  localized  region  near  the  impact  site,  and  could  be  incorporated  into 
the  contact  model  rather  than  the  stiffness  of  the  laminate.  This  allows 
the  initial  stiffness  of  the  laminate  to  be  used  throughout  the  analysis. 
This  assumption  is  validated  to  some  extent  by  experiments  performed  by  AFML 
on  plates  with  either  epoxy  or  thermoplastic  matrices. 

Figures  1  and  2  show  the  experimental  responses  of  graph ite/epoxy  and 
graphite/PEEK  plates  to  the  same  impact  event.  Figure  1  shows  the  contact 
force  history  for  8  ply  laminates  simply  supported  over  a  50  mm  diameter 
circle,  and  Figure  2  shows  the  same  for  16  ply  laminates.  The  two  matrix 
systems  used  in  these  plates  have  essentially  the  same  initial  stiffness, 
but  the  thermoplastic  is  much  more  ductile.  The  figures  show  that  the 
response  of  the  two  material  systems  are  similar  up  to  the  initiation  of 
damage,  which  is  consistent  with  the  assumption  stated  above.  It  is  further 
assumed  that  differences  in  the  point  at  which  damage  occurs  may  be 
attributed  to  differences  in  the  strength  of  the  two  matrices,  rather  than 
differences  in  the  dynamic  response. 

Objectives  of  the  impact  analysis  are  to  predict  both  an  accurate 
flexural  response  and  an  accurate  stress  distribution.  The  level  of 
refinement  needed  to  predict  accurate  stresses  in  a  finite  element  model  is 
typically  much  greater  than  the  refinement  necessary  for  accurate 
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displacements.  This  results  in  a  dynamic  model  which  includes  for  more 
degrees  of  freedom  than  are  necessary  for  predicting  displacement  histories. 

To  overcome  this  problem,  the  laminate  is  modeled  as  a  number  of 
substructures,  each  of  which  is  representative  of  a  smaller  portion  of  the 
complete  laminate.  A  reduction  procedure  is  used  which  lefines  the 
effective  mass  and  stiffness  of  the  substructure  in  terms  of  a  small  number 
of  "master"  degrees  of  freedom.  The  response  of  the  complete  substructure 
may  then  be  approximated,  given  the  response  of  the  master  degrees  of 
freedom. 

The  reduced  substructures  are  then  assembled,  similar  to  the  assembly 
of  typical  finite  elements,  into  a  model  of  the  complete  laminate.  The  same 
type  of  reduction  is  then  performed  on  the  full  model,  allowing  the  response 
of  the  complete  set  of  master  degrees  of  freedom  to  be  approximated  by  the 
response  of  a  set  of  "governor"  degrees  of  freedom. 

Adding  the  contact  model  to  the  reduced  dynamic  model  of  the  laminate 
then  gives  the  reduced  equations  of  motion  for  the  plate/impactor  system. 
These  equations  are  integrated  to  give  the  response  of  the  governor  degrees 
of  freedom.  This  response  is  expanded  to  give  the  response  of  all  master 
degrees  of  freedom,  which  are  then  used  to  give  the  complete  displacement 
response  of  each  substructure.  This  model  building  sequence  is  shown 
schematically  in  Figure  3. 

For  a  given  substructure,  the  displacement  solution  at  a  point  in  time 
is  used  to  calculate  strains.  From  this  strain  distribution  stresses  are 
calculated  and  used  in  a  suitable  failure  criteria  to  assess  failure  in  that 
substructure  at  that  time.  These  calculations  are  performed  for  each 
substructure  at  regular  time  intervals  to  predict  the  onset  of  damage. 

The  analytical  procedure  then  consists  of  three  phases:  substructure 
generation  and  reduction;  global  model  assembly,  reduction,  and  solution; 
and  back  substitution,  stress  and  failure  analysis.  Since  the  mass  and 
stiffness  of  the  laminate  are  assumed  to  be  constant,  the  model  generation, 
assembly  and  reduction  can  be  performed  once  and  used  throughout  the 
analysis . 
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2.  LAMINATED  SOLID  FINITE  ELEMENTS 


A  number  of  displacement  shape  functions  were  considered  for  use  in 
this  program.  Two-dimensional  elements  were  developed  with  these  shape 
functions  and  used  to  evaluate  the  relative  merits  of  the  functions.  An 
element's  ability  to  approximate  the  complex  stress  distributions  through 
the  thickness  of  a  short  beam  in  three  point  bending  was  used  as  a 
benchmark. 

From  these  analyses,  it  was  determined  that  the  most  appropriate 
element  for  further  development  was  a  subparametric  element  which  utilized 
cubic  displacement  shape  functions,  and  included  first  derivatives  of 
displacements  as  nodal  variables.  The  details  and  findings  of  this  study 
have  been  documented  in  Reference  2. 

Cubic  shape  functions  are  particularly  well  suited  for  modeling  beams 
and  plates  since  the  classical  beam  and  plate  theories  assume  cubic 
displacement  fields.  Shear  strain  distributions  in  this  element  are  then 
quadratic,  allowing  a  parabolic  distribution  to  be  modeled  with  a  single 
element  through  the  thickness.  More  complex  distributions  are  approximated 
with  only  a  few  elements  through  the  thickness. 

The  inclusion  of  derivatives  as  nodal  variables  is  also  attractive  in 
that  it  allows  boundary  conditions  to  be  defined  In  terms  of  strain.  For 
example,  these  derivatives  may  be  used  in  conjunction  with  the  constitutive 
law  for  a  surface  ply  to  inforce  stress  boundary  .mditions  at  the  top  and 
bottom  surfaces  of  a  laminate. 

A  two-dimensional  element  was  developed  for  modeling  beams,  and  three- 
dimensional  elements  were  developed  for  modeling  plates.  Both  types  of 
elements  were  incorporated  into  computer  programs  which  generate  and  reduce 
beam  or  plate  substructures,  as  previously  described.  Descriptions  of  these 
programs  are  included  in  the  Appendix. 

2-D  Element 

The  two-dimensional  element  is  based  upon  the  work  of  Webster  (3), 
extended  to  account  for  step  changes  in  material  properties  through  the 
thickness  of  the  element.  This  extension  allows  a  single  element  to  model 
several  plies  in  a  laminate. 
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A  subparametric  fonaulation  is  used,  In  which  the  element  geometry  is 
restricted  to  be  a  unit  width  rectangle,  with  sides  parallel  to  the  global 
axes.  The  2*D  geometry  and  nomenclature  are  shown  in  Figure  4.  Geometric 
shape  functions  for  the  elements  are  given  by 

4 

g(x.z)  -  Z  N.g.  (1) 

i-1  1  1 


where 

gi  -  g(xi,zi) 

»!  *  t  <l«0>(1+2o> 

and 

X  -  XX . 
o  i 


where  (x^,z^)  are  the  local  coordinates  of  node  i. 

Nodal  variables  for  the  element  consist  of  translations  in  the  x  and  z 
directions  (u  and  w)  ,  and  the  first  derivatives  of  these  displacements  with 
respect  to  x  and  z  (u,x,  u>z,  w,x,  and  w>z).  The  displacement  shape 
function  is 

4 

f (x , z)  -  Z  [N]  If  )  (2) 

i-1 


where 

T 

lf>I- 

tfi 

f.  f, 

Xi  Zi 

[N]i  - 

(Ni 

N  N  ] 

Xi  Zi 

and  the  components  of  [N]^  are  given  by 
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(3a) 


Ni  "  8(1+xo)(1+Zo)(2+Xo(1*Xo)  +  zo(1'2o>) 

Nx  “  *8(1+xo)2(1+Zo)(1'Xo)xi  (3b) 

i 

N  -  -£(l+x  >(l+z  )2(l-z  )z.  (3c) 

Z .  O  O  O  O  1 

1 

Strain  components  which  contribute  to  the  strain  energy  in  the  element 

are  the  extensional  strains  t  and  t  ,  and  the  shear  strain  7  .  These 

X  z  zx 

strains  are  given  in  terms  of  the  displacements  as 


Using  the  shape  function  given  by  (2),  the  strain  in  the  element  may  be 
expressed  in  terms  of  the  nodal  variables 

U 

{ t }  -  Z  [B].  (u)  -  [B]{u)  (5) 

i-1  1  x 

where 


6 


The  stress  in  the  element  is  given  by 


where  [C]  is  the  constitutive  law  that  has  been  modified  to  reflect  a  state 
of  either  plane  stress  or  plane  strain.  Using  (5),  the  stress  may  be 
expressed  in  terms  of  the  nodal  variables  as 

[a)  -  [C] [B] (u)  (7) 

Using  (5)  and  (7),  the  stiffness  of  the  element  may  be  found  by 
minimizing  potential  energy,  yielding 

[Ke]  -  JV[B]T[C][B]dv  (8) 

Since  the  material  properties  in  a  laminate  are  piecewise  constant  through 
the  thickness,  the  integral  in  (8)  may  be  expressed  as 

n  i/2  z\  _ 

[K  1  -  E  f  f  [B] A [C]  [B]  dzdx  (9) 

i-1  -i/2  zj  1 

where  z*  and  z^  are  the  locations  of  the  top  and  bottom  surfaces  of  the  ith 
ply  respectively,  and  [C]^  is  the  constitutive  law  for  the  ith  ply. 

A  consistant  mass  matrix  is  generated  for  the  element,  using  the 

geometric  shape  function  given  by  (1).  This  matrix  is  evaluated  using  the 
Integral 

-  fV  (N']Tp[N']dv  (10) 

As  with  the  constitutive  law,  the  density  of  the  laminate  is  piecewise 
constant  through  the  thickness.  The  mass  matrix  may  then  be  evaluated  by 
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(11) 


n  i/2  zt  T 

[Ml  -  T  J  J  [N']  V(N']dzdx 

e  i-1  -i/2  z*  1 

where  p ^  is  the  density  of  the  ith  ply. 

This  mass  matrix  only  associates  inertia  with  translational  degrees  of 
freedom.  However,  since  there  is  a  mass  distribution  through  the  thickness, 
an  inertial  effect  analogous  to  rotary  inertia  is  included. 

When  a  displacement  solution  has  been  found,  the  strain  at  a  point  in 
the  element  may  be  evaluated  using  (5).  The  stress  corresponding  to  this 
strain  is  then  given  by  (6).  Average  strain  in  a  ply  may  be  found  by 
integrating  (5)  over  the  volume  of  that  ply  contained  in  the  element.  This 
average  strain  is  given  by 

l  i/2  2t 

{«}.  “  /  J  [B] dzdx) (u)  (12) 

i  i/2  z[ 


where  is  the  volume  of  ply  i  contained  in  the  element.  The  average 
stress  in  the  ply  is  then  given  by 

l^li  -  [C]i{i}i  (13) 

A  similar  procedure  may  be  performed  to  determine  the  average  states  of 
stress  and  strain  at  a  ply  interface.  It  is  assumed  that  the  interface  has 
no  thickness,  allowing  the  Integration  to  be  performed  over  the  interface 
area,  rather  than  volume.  Using  (5),  the  average  strain  at  an  interface  is 


1  i/2 

(T*  J  (B(x,z.) ]dx) (u) 
j  -i/2  J 


(14) 


where  A^  is  the  surface  area  of  the  jth  interface, 
The  average  stress  at  the  interface  is  then 


contained  in  the  element. 


(ojj  -  [CJjIOj 


(15) 
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where  [C)j  is  the  constitutive  law  for  the  interface.  Note  that  since  the 
volume  of  the  interface  is  zero,  it  does  not  contribute  to  the  mass  or 
stiffness  of  the  element. 


3-E  Eltawtg 

Two  three-dimensional  elements  were  developed,  both  of  which  utilize 
cubic  displacement  functions  and  retain  first  derivatives  as  nodal  degrees 
of  freedom.  The  first  of  these  is  an  extension  of  the  two-dimensional 
element,  and  requires  that  the  element  geometry  be  a  rectangular  hexadedron 
with  sides  parallel  to  the  global  coordinate  axes  (Figure  5).  The  second 
element  is  less  restrictive  in  that  it  allows  the  geometry  to  be  a  uniform 
thickness  prism,  defined  by  an  arbitrary  quadrilateral  in  the  x-y  plane 
(Figure  6) .  Since  the  rectangular  element  is  a  special  case  of  the 
prismatic  element,  the  formulation  of  the  general  geometry  will  be  given 
first,  and  degenerated  for  the  special  case. 

Formulation  of  the  element  is  simplified  by  mapping  the  element 
geometry  onto  a  two  unit  cube  as  shown  in  Figure  7.  Edges  of  the  element 
are  restricted  to  be  straight  lines,  allowing  a  subparametric  formulation. 
In  the  local  (r,s,q)  coordinate  system,  the  geometric  shape  function  is 
given  by 

8 

g(r,s,q)  -I  Nj  g.  (16) 

i-1  1  1 


where 


gi  "  g(ri*  si’  qi) 

”1  -  8<1+ro><UVa+V 

and 


r 

o 


rr 


i 


s  -  ss 
o 


q0  -  qq* 


where  (r^,s^,q^)  are  the  coordinates  of  node  i  in  the  local  coordinate 
system. 

Nodal  variables  for  this  element  consist  of  translations  in  each 
coordinate  direction  (u,v  and  w) ,  and  the  first  derivatives  of  these 
translations  with  respect  to  each  coordinate.  Displacement  shape  functions 
are  cubic  polynomials  in  the  local  coordinate  system,  given  by 

2  2  2 

f(r,s,t)  -  cq  +  CjT  +  c2s  +  c3q  +  c4r  +  c5rs  +  cgs  +  c?sq  +  cgq  (17) 

3  2  2  3  2  2 

+  cgqr  +  c1Qr  +  cur  s  +  c^rs  +  c13s  +  c14s  q  +  c^sq 

3  2  2  3  3  3 

+  c16q  +  c1?q  r  +  c^qr  +  c^rsq  +  c2Qr  s  +  c21rs  +  c^s  q 

3  3  3  2  2  2 

+  c23sq  +  c24q  r  +  c^qr  +  c^r  sq  +  c^rs  q  +  c^rsq 

3  3  3 

+  c29r  sq  +  c3Qrs  q  +  c31rsq 

Expressed  in  terms  of  the  local  nodal  variables,  the  displacement  shape 
functions  are 
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f(r.s.q)  -  S  [N] . If  ) 
i-1  1 


where 


and 


(18) 


INlf  -  (Nt  N 


N 


N  ] 

Si  V 


where  the  components  of  [NJ^  are  given  by 


»,  - 
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(1+r  )(l+s  )(l+q  )(2+r  (1-r  )  +  s(l-s)  +  qo(l-qn>) 


"rj  -  -U  <1«„>2<lt‘o)<l+’o><1-ro>rl 


(19a) 

(19b) 
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(19c) 


\  -  *16  a«0)(i«.)2(1+%)<l'so)si 


16 


(196) 


The  nodal  derivatives  used  in  (18)  are  taken  with  respect  to  the  local 
coordinates,  and  oust  be  transformed  into  the  global  coordinate  system. 

Using  the  geometrical  shape  functions  given  in  (16),  and  the  chain  rule,  the 
transformation  from  global  to  local  variables  is  given  by 


(f ' )±  “  [J'litfli 
where 


(20) 


-  (f± 
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Substituting  (20)  into  (18)  yields  the  displacement  shape  function  in 
terms  of  the  global  nodal  variables 


0  o 

f(r,s,q)  -  E  [N] . [ J  ]  (f).  -  E  (N*].{f). 

i-1  1  i-1  11 


(21) 
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A  complete  three-dimensional  state  of  strain  is  assumed  for  the  element 


(€}--  cz  >  "  dz  (22) 

fiy  £u 

7xy  3x  3y 

£w  3v 

7yz  3y  32 

au  +  fiw 

7zx  3z  3x  v 

The  partial  derivatives  in  (22)  are  evaluated  using  (21),  to  give 
strains  in  terms  of  the  nodal  variables.  Formulation  of  the  element  mass 
and  stiffness  matrices  follow  the  same  procedures  used  for  the  two- 
dimensional  element.  Similarly,  average  stresses  and  strains  may  be  found 
by  integrating  over  the  volume  of  a  ply,  or  the  area  of  an  interface. 

By  restricting  the  geometry  of  the  element  to  be  rectangular, 
formulation  of  the  element  matrices  can  be  greatly  simplified.  Shape 
functions  for  this  element  may  be  found  by  substituting  the  identities 

r  -  2 x/i 

s  -  2y/w  (23) 

q  -  2z/t 

into  (16)  and  (19).  The  transformation  from  local  to  global  nodal  variables 
then  becomes 

f£i  .  i  f£i 

3r  2  3x 
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l£i  _tiii 

at  2  dz 

The  computational  effort  required  for  this  special  case  is 
significantly  reduced  relative  to  that  required  for  the  general  element.  As 
with  the  two-dimensional  element,  consistent  mass  matrices  are  generated  for 
both  three-dimensional  elements. 

3.  NON-LINEAR  CONTACT  MODEL 

The  impactor  is  modeled  as  a  lumped  mass,  coupled  to  the  laminate 
through  a  nonlinear  spring  at  a  single  node.  This  model  is  based  on  the 
work  of  Sun  &  Yang,  and  Yang  &  Sun  (Ref.  4  and  5).  The  contact  force  is 
assumed  to  be  a  nonlinear  function  of  the  difference  in  the  lateral 
displacements  of  the  impactor  and  the  impacted  node. 

When  loading,  the  contact  force  is  given  by 

qx  -  KAn  (24) 

when  A  is  the  difference  in  the  lateral  deflections,  and  K  and  n  are 
empirically  defined  loading  parameters. 

When  unloading,  the  contact  force  is  given  by 


*>1  ■  V 


A  -  A 


A  -  A 
max  o 


where  m  is  an  empirical  parameter,  and  Aq  represents  a  permanent  indentation 
in  the  laminate  given  by 


A  -  A  [1  -  (t“  )  ] 
o  max1  A  J 
max 


(26a) 


when  A  is  greater  than  A  ,  and  by 
max  &  cr’  J 


Ao-0 


(26b) 
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which  are  now  the  reduced  equations  of  motion  for  the  complete  laminate. 

The  effect  of  the  reduction  process  is  to  lump  all  inertial  effects  at  the 
governor  degrees  of  freedom.  Great  care  should  then  be  exercised  in 
selecting  both  master  and  governor  degrees  of  freedom. 

The  equations  of  motion  for  the  complete  plate/impactor  system  requires 
the  equation  of  motion  for  the  impactor.  Since  the  iapactor  is  included  as 
a  lumped  mass,  its  equation  of  motion  is  given  by  rigid  body  dynamics  as 

“luI  "  qI  (33) 

where  m^  is  the  mass  of  the  impactor,  u^  is  its  displacement,  and  qj  is  the 
contact  force. 

The  set  of  simultaneous  differential  equations  given  by  (32)  and  (33) 
may  be  integrated  to  give  the  dynamic  response  of  the  complete  system.  This 
integration  is  most  easily  performed  by  transforming  the  second  order 
differential  equations  into  two  sets  of  first  order  equations.  Defining 

(ug)  -  (vg)  (34a) 

and 

Uj  -  Vj  (34b) 

Equations  (33)  and  (32)  may  be  solved  for  (ug)  and  u^,  yielding 
Vj  -  ui  -  qj/Oj  (35a) 

and 

(v  )  -  <V  “  lM*l  'X(tq*>  *  [ K* ]  { u  } )  (35b) 

o  ©  © 


Combining  (34a),  (34b) 


(35a),  and  (35b)  gives 


{u' }  -  ( q ' )  +  [D' ] (u' ) 


(36) 


where 

.  -  .T  ,  T  T, 

{U  )  -  lUj  Ug  Vj  vg) 

[D']  -  0  0  1  0 

0  0  0  [I] 

0  0  0  0 

0  - [M*] * 1 [K*]  0  0 

and 

( q ”  )T  -  (0  (0)T  qj/nij  ([M*]'1{q*))T) 

Equation  (36)  can  easily  be  integrated  numerically  using  one  of  the 
Runge -Kutta  algorithms,  which  'step  out’  the  dynamic  response  of  the  system, 
given  initial  displacements  and  velocities.  The  general  form  of  these 
algorithms  for  the  solution  at  the  nth  time  step  is 

{u')n  -  {u ')nl  +  At*  f ( t ,u  ,  u  )  (37) 

^  9  f 

where  At  is  the  time  increment,  and  f(t,u  ,u  )  is  an  estimate  of  the  average 
slope  over  the  time  interval.  Explicit  forms  of  Equation  (37)  are  given  in 
Reference  7. 

Using  (37),  the  displacements  and  velocities  of  the  impactor  and 
governor  degrees  may  be  stepped  out  to  any  point  in  time.  Accelerations  at 
this  time  are  given  by  (35a)  and  (35b).  The  response  of  the  complete  set  of 
master  degrees  of  freedom  are  given  by  (31),  allowing  the  response  of  the 
slave  degrees  of  freedom  to  be  found  using  (28). 
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FAILURE  ANALYSIS 


The  displacement  history  given  by  (37),  (31),  and  (28)  may  be  used  in 
Equations  (12)  and  (14)  to  determine  average  ply  and  interface  strains  in 
each  element.  These  strains  are  then  used  in  Equation  (13)  and  (15)  to 
determine  the  corresponding  stresses.  Thus,  the  time  varying  spatial 
distributions  of  stress  and  strain  are  given  throughout  the  laminate. 

Damage  initiation  may  then  be  determined  by  applying  failure  criteria 
to  either  the  stress  or  strain  distributions.  The  first  criterion  that  is 
meet  will  determine  the  time  and  location  of  the  onset  of  damage.  The 
predicted  response  after  the  initiation  on  damage  is  then  of  questionable 
value  since  it  will  not  reflect  changes  in  stiffness  which  may  effect  the 
actual  response.  This  will  be  discussed  in  further  detail  in  the  analysis 
section  of  this  report. 

Since  all  six  stress  components  are  available,  it  is  desirable  to  use 
failure  criteria  that  account  for  interactions  between  these  components. 
For  predicting  failure  within  a  ply,  Hashin's  failure  criteria  with  slight 
modifications  are  used  (Ref.  8).  These  criteria  define  distinct  modes  of 
failure  for  both  the  fiber  and  matrix. 

The  fiber  tensile  mode  assumes  a  quadratic  interaction  between  the 
axial  stress  and  the  maximum  axial  shear  and  is  predicted  when 


,  2 
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12 
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31 
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(38) 


where  s^  and  s^  are  the  axial  tensile  and  axial  shear  strengths 
respectively.  This  criterion  is  applicable  when  is  positive,  and  is  used 
to  predict  a  failure  mode  characterized  by  fiber  breakage. 

Vhen  a ^  is  compressive  it  is  assumed  that  failure  is  characterized  by 
fiber  buckling  and  is  only  dependent  upon  The  compressive  fiber  mode 

failure  criterion  is  then  given  by  the  maximum  stress  criterion 


(39) 


where  s^  is  the  axial  compressive  strength  of  the  pi-'. 
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Matrix  node  failures  are  characterized  by  cracks  running  parallel  to 
the  fibers.  Failures  are  described  as  being  tensile  or  compressive, 
depending  upon  the  sign  of  the  quantity  ( o ^  +  o^).  Both  natrix  mode  failure 
criteria  assume  quadratic  interactions  between  the  transverse  stresses,  both 
in-plane  and  through  the  thickness,  the  maximum  shear  in  the  transverse 
plane  and  the  maximum  axial  shear. 

When  (a 2  +  o^)  is  positive,  the  tensile  mode  criterion  is  used.  This 
criterion  is  given  by 


°2  *  °3 


(t?3  *  ‘W  ^  (fl?  *  ill? 

n  '  a 


-  1 


(40) 
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where  s£  and  S2 3  are  the  transverse  tensile  and  transverse  shear  strengths. 
For  (oj  +  ^3)  negative,  the  compressive  failure  criteria  is  given  by 


+  r 


’23 


2  2  2 
23  (r 12  +  r31) 

^  r%  "  ^ 
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(41) 


which  is  a  simple  quadratic  interaction  between  the  maximum  transverse  and 
axial  shears.  Failure  predicted  by  this  criterion  will  be  referred  to  as 
shear  failures. 

In  contrast  to  the  plies,  interfaces  are  assumed  to  be  homogeneous, 
matrix  rich  regions.  Although  orthotropic  interface  properties  could  be 
included  in  the  problem  formulation,  commonly  used  matrix  materials  are 
typically  isotropic,  allowing  the  use  of  simpler  failure  criteria. 

Assuming  an  isotropic  interface,  a  simple  criterion  may  be  established 
by  postulating  a  failure  surface  that  is  analogous  to  the  Von  Mises  yield 
surface.  Failure  is  then  predicted  by  a  maximum  stress  criterion  which 
utilizes  the  octahedral  shear  stress  (Ref.  9) 

\^ala2>2  +  (<72*<73>2  +  <ct3'<7i)2]/t  “  1  (42) 


where  a ^  are  principle  stresses,  and  r  is  the  shear  strength  of  the 
interface . 
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It  should  be  noted  that  a  number  of  criteria  have  been  postulated  for 
predicting  the  failure  of  fibrous  composites.  The  set  of  criteria  described 
above  constitute  only  one  of  these.  Its  applicability  and  the  applicability 
of  any  such  criteria  depend  upon  the  characteristics  of  the  materials  to 
which  they  are  applied,  and  can  only  be  validated  experimentally. 

This  set  of  criteria  is  then  presented  as  one  possible  set  which  has 
been  shown  to  agree  well  with  experimental  results  for  epoxy  composites. 
Other  criteria  may  be  better  suited  for  different  materials  and  could  be 
substituted  into  the  framework  of  this  analysis. 
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III.  ANALYSIS 


A  number  of  laminated  composite  plates  were  analyzed  to  verify  the 
applicability  of  the  analytical  model.  These  analyses  included  AS4/3502 
plates  that  were  tested  by  AFML,  and  IM6/934  plates  that  were  part  of 
another  research  program  being  performed  by  MSC.  The  AFML  data  were 
collected  for  a  number  of  different  support  geometries  and  laminate 
configurations,  whereas  the  IM6/934  analyses  considered  identical  plates 
impacted  at  different  velocities.  The  experimental  data  consisted  of 
contact  force  histories  for  the  impact  events. 

1 .  MODEL  DESCRIPTIONS 

Experimental  impact  data  for  quasi -isotropic  AS4/3502  plates  containing 
either  8  or  16  plies  was  made  available  by  AFML.  These  plates  were  simply 
supported  over  25,  50,  75  or  100  mm  circles,  and  a  pendulum  impactor  was 
used  for  all  tests.  This  impact  test  setup  is  described  in  Reference  10. 

Of  these  plates,  the  50  and  75  mm  support  16  ply,  and  the  25,  50  and  75  mm 
support  8  ply  were  analyzed.  Table  1  outlines  the  conditions  for  each  of 
these  analyses. 

Ply  properties  for  AS4/3502  were  taken  from  References  11  and  12,  and 
neat  resin  properties  for  3502  were  taken  from  13.  These  properties  are 
summarized  in  Table  2  and  were  used  for  all  AS4/3502  analyses. 

Parameters  for  the  contact  model  were  taken  from  experimental  tests  on 
a  graphite/epoxy  laminate  conducted  by  Sun  and  Wang  (14).  Although  the 
laminate  used  for  these  tests  was  not  identical  to  the  AS4/3502  laminate 
being  considered  here,  it  was  assumed  that  the  contact  parameters  for  the 
two  would  be  similar.  The  parameters  used  for  the  AS4/3502  analyses  are 
summarized  in  Table  3. 

Finite  element  models  of  both  the  8  and  16  ply  laminates  were  generated 
using  the  3-D  laminated  solid  element  which  allows  arbitrary  geometries  in 
the  X-Y  plane.  These  models  were  constructed  such  that  they  could  be 
constrained  over  a  25,  50  or  75  mm  circle,  as  shown  by  the  X-Y  mesh  for  a 
quarter  of  the  laminate  in  Figure  9.  Both  models  utilized  this  mesh,  and 
included  two  elements  through- the -thickness. 
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Two  quarter-plate  substructures  were  generated,  reduced  and  assembled 
into  a  dynamic  model  of  half  the  plate.  Displacement  boundary  conditions 
were  then  added  to  this  model  to  simulate  the  circular  support,  and  to 
enforce  rotational  symmetry  about  the  center  of  the  plate.  The  reduced 
dynamic  model  included  39  governor  degrees  of  freedom. 

The  IM6/934  plates  were  quasi-isotropic  32  ply  laminates,  simply 
supported  over  a  76.2  mm  X  127.  mm  rectangle  as  summarized  in  Table  4. 

These  tests  utilized  a  verticle  drop  tower,  in  which  the  height  of  the 
inpactor  was  varied. 

The  finite  element  model  generated  for  these  analyses  utilized  the  3-D 
laminated  solid  element  that  is  restricted  to  rectangular  geometries. 
Development  of  the  dynamic  model  for  these  analysis  was  similar  to  that  for 
the  AS4/3502  model  where  the  rotational  symmetry  of  the  laminate  was 
exploited.  The  quarter -plate  mesh  for  the  IM6/934  model  is  shown  in  Figure 
10. 

2 .  RESULTS 

The  first  AS4/3205  plate  analyzed  was  a  16  ply  laminate  supported  over 
a  50  mm  diameter  circle.  The  response  of  this  plate  was  stepped  out  over 
the  first  millisecond  of  the  impact  event  using  a  time  step  of  50 
nanoseconds.  The  displacement  solution  was  expanded,  and  a  stress  analysis 
was  performed  every  20  microseconds. 

Lateral  displacement  histories  predicted  for  the  impactor  and  impacted 
node  are  shown  in  Figure  11.  As  indicated  by  the  slope  of  the  impactor 
displacement  history  in  this  figure,  the  velocity  of  the  impactor  is  only 
slightly  reduced  over  this  interval.  Very  little  of  the  kinetic  energy  of 
the  impactor  has  been  transferred  into  the  plate  at  this  point. 

The  difference  between  the  two  displacement  histories  in  Figure  11 
defines  the  magnitude  of  the  contact  force.  Al thought  it  is  not  easily  seen 
from  the  figure,  there  is  a  high  frequency  displacement  component  contained 
in  the  history  for  the  impacted  node.  This  component  is  more  readily  seen 
in  the  contact  force  history  which  is  compared  to  the  experimental  contact 
force  in  Figure  12. 

Good  agreement  is  seen  between  the  experimental  and  analytical  contact 
force  histories.  Of  particular  interest  is  the  aforementioned  high 
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frequency  component  that  has  approximately  the  same  frequency  and  peak  to 
peak  magnitude  in  both  cases.  This  correlation  lends  credibility  to  the 
model  as  a  whole,  and  justifies  the  use  of  the  chosen  contact  parameters. 

It  is  argued  that  the  abrupt  drop  off  in  the  experimental  contact  force 
at  about  0.8  ms  indicates  the  initiation  of  damage  in  the  laminate.  This 
point  in  the  contact  history  then  serves  as  a  reference  for  assessing  the 
models  ability  to  predict  the  onset  of  damage. 

The  stress  histories  predicted  by  this  analysis  indicate  tensile  matrix 
failures  developing  as  soon  as  0.4  ms  into  the  event.  This  damage  is  found 
in  the  bottom  most  ply  of  elements  closest  to  the  impactor.  It  was  assumed 
that  these  failures  would  result  in  cracks  in  the  matrix  that  would  not 
significantly  effect  the  stiffness  of  the  laminate,  and  the  analysis  was 
continued.  The  validity  of  this  assumption  will  be  discussed  in  greater 
detail  later  in  this  report. 

As  the  analysis  progresses,  the  region  experiencing  tensile  matrix 
failure  increases.  This  continues  until  about  0.92  ms  when  a  compressive 
matrix  mode  (shear)  failure  is  predicted  for  the  top  ply  in  the  elements 
nearest  the  impactor.  From  this  point  on,  these  shear  failures  quickly 
progress  down  toward  the  center  of  the  laminate  as  tensile  failures  progress 
up.  The  validity  of  the  solution  after  0.92  ms  is  then  very  doubtful. 

The  development  of  the  stress  state  which  produces  the  shear  failure  is 
shown  in  Figure  13  which  presents  histories  for  the  maximum  transverse  and 
axial  shears  in  the  top  ply.  From  this  figure  it  can  be  seen  that  although 
the  failure  is  dominated  by  the  maximum  transverse  shear,  there  is  still  a 
substantial  contribution  from  the  axial  shear.  If  it  is  assumed  that  this 
failure  produces  the  damage  evident  in  the  experimental  contact  force  trace, 
the  difference  in  the  predicted  to  experimental  damage  initiation  times  is 
on  the  order  of  15%.  Given  the  uncertainty  in  shear  strengths,  this  is 
fairly  good  agreement. 

If  this  damage  initiation  scenario  is  accepted,  it  becomes  necessary  to 
refine  the  definition  of  damage  initiation.  That  is,  damage  initiation  may 
be  taken  as  the  first  failure  that  is  detectable  upon  examination  of  the 
laminate,  or  the  first  failure  that  is  detectable  through  the  response  of 
the  laminate.  These  two  failures  may  or  may  not  be  one  and  the  same. 

Examining  the  IM6/934  analysis  helps  clarify  this  point.  In  this 
analysis,  a  quasi  -  isotropic  plate  was  impacted  at  a  velocity  high  enough  to 
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develop  tensile  cracks,  but  below  that  necessary  for  the  initiation  of  shear 
failures.  The  analytical  contact  force  history  predicted  by  this  analysis 
is  shown  in  Figure  14,  along  with  the  point  at  which  damage  initiation  was 
predicted. 

As  with  the  AS4/3502  plate,  this  initial  damage  is  a  tensile  matrix 
mode  failure  in  the  bottom  most  ply.  Proprietary  experimental  results  for 
this  plate  showed  no  evidence  of  damage  in  the  contact  force  history.  Post 
examination  of  the  plate  did,  however,  show  visable  cracks  on  the  back 
surface.  Damage  was  then  initiated  in  the  plate,  but  did  not  effect  its 
dynamic  response. 

This  plate  was  also  analyzed  at  higher  impact  velocities,  for  which  the 
damage  initiation  scenario  was  similar  to  that  of  the  IM6/3502  plate.  These 
analyses  correlated  very  well  with  proprietary  experimental  results  for  both 
the  contact  force  history  and  the  damage  initiation  times. 

These  analyses  both  verified  the  existence  of  tensile  cracks  early  in 
the  dynamic  response,  and  validated  the  use  of  the  analytical  model  after 
these  cracks  were  predicted.  Analysis  of  the  AS4/3502  plates  was  then 
resumed. 

The  next  case  considered  was  a  16  ply  plate  supported  over  a  75  mm 
circle.  The  response  of  this  plate  was  similar  to  the  50  mm  support  case, 
and  followed  the  same  damage  scenario.  The  experimental  and  analytical 
contact  force  histories  for  this  case  are  compared  in  Figure  15.  As  shown 
in  the  figure,  damage  consists  of  the  formation  of  tensile  cracks  followed 
by  shear  failures  in  the  top  ply.  Throughout  the  remainder  of  this  report 
damage  initiation  will  be  interpreted  as  the  initiation  of  failures  which 
influence  the  dynamic  response. 

Eight  ply  AS4/3502  plates,  supported  over  25,  50  and  75  mm  circles, 
were  the  final  analyses  performed.  Experimental  and  analytical  contact 
force  histories  for  these  cases  are  shown  in  Figures  16,  17  and  18 
respectively.  In  all  cases,  the  formulation  of  damage  was  similar  to  that 
seen  in  the  16  ply  plates. 

Analytical  results  for  the  25  mm  support  case  correlate  well  with  the 
magnitude  of  the  experimental  contact  force  up  to  the  time  of  damage 
initiation.  The  50  and  75  mm  support  cases,  however,  only  correlate  through 
the  first  millisecond  of  the  response  and  then  diverge.  This  implies  that 
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there  Is  considerable  membrane  stretching  in  these  plates,  a  non-linear 
effect  that  is  not  included  in  the  formulation  of  the  model. 

This  is  verified  by  examining  the  magnitudes  of  the  extensional  strains 
predicted  by  the  linear  formulation.  Using  the  displacement  solution  for 
the  node  on  the  back  surface  of  the  plate  at  X-  2.54  am,  Y  -  0.0  mm,  global 
strains  were  calculated  using  both  the  linear  and  non-linear  definitions  of 
strain.  Since  first  derivatives  of  the  displacements  are  included  as  nodal 
degrees  of  freedom,  the  linear  strains  are  found  directly,  and  the  non¬ 
linear  strains  are  easily  calculated. 

The  dominant  extensional  strain  at  this  point  is  the  strain  in  the  Y 
direction.  Histories  for  the  linear  and  non-linear  formulations  of  this 
component  are  shown  in  Figure  19.  As  expected,  the  two  curves  do  diverge. 

It  should  be  noted,  however,  that  the  non-linear  curve  is  calculated  from  a 
displacement  solution  predicted  by  the  linear  theory.  A  different  curve 
will  result  if  the  non-linear  theory  is  used  throughout  the  analysis,  but  a 
similar  conclusion  is  expected. 
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The  analyses  included  in  the  program  have  helped  to  define  the  limits 
of  applicability  of  the  analytical  model,  and  have  shown  it  to  be  an 
accurate  engineering  tool  when  used  within  these  limits.  The  model  is 
particularly  well  suited  for  the  analysis  of  plates  with  low  aspect  ratios 
which  typically  exhibit  significant  shear  deformation.  This  deformation  can 
be  characterized  very  well  by  the  higher  order  (cubic)  displacement  shape 
functions  used  in  the  laminated  solid  elements. 

Thinner  plates  will  exhibit  stiffening  due  to  membrane  stretching,  much 
like  the  stiffening  of  a  drum  head  as  it  is  tightened.  This  non-linear 
effect  can  only  be  included  by  utilizing  a  higher  order  definition  of  strain 
as  opposed  to  a  higher  order  shape  function.  Although  the  use  of  the  non¬ 
linear  strain  tensor  in  a  finite  element  is  well  defined,  in  practice,  it 
presents  severe  problems  since  it  necessitates  reformulating  the  stiffness 
matrix  periodically.  If,  however,  thin  plates  are  to  be  analyzed,  the  non¬ 
linear  formulation  must  be  considered. 

By  qualitatively  comparing  the  aspect  ratio  of  each  plate  analyzed  to 
the  degree  of  correlation  between  the  experimental  and  analytical  contact 
force  histories,  an  upper  limit  on  the  aspect  ratio  can  be  determined. 

Table  5  lists  the  aspect  ratios  for  each  of  the  plates  considered  in  this 
program.  For  circular  plates,  the  aspect  ratio  is  taken  as  the  ratio  of  the 
diameter  to  the  thickness;  for  the  rectangular  plate,  major  and  minor  aspect 
ratios  are  given  by  the  ratios  of  the  major  and  minor  support  lengths  to  the 
thickness  respectively. 

Good  agreement  was  found  between  the  experimental  and  analytical 
contact  force  histories  for  the  50  and  75  mm  support,  16  ply  AS4/3502 
plates;  the  25  mm  support,  8  ply  AS4/3502  plate;  and  the  IM6/934  plate.  All 
of  these  plates  had  aspect  ratios  less  than  40.  The  remaining  plates  had 
aspect  ratios  greater  than  40,  and  did  not  show  as  good  an  agreement  with 
the  experimental  data.  An  aspect  ratio  of  about  40  is  then  assumed  to  be 
the  upper  limit  for  which  the  model  is  appropriate. 

All  analyses  predicted  matrix  cracking  in  plies  on  the  tensile  side  of 
the  laminate  early  in  the  response.  The  existence  of  these  cracks  has  been 
verified  in  unpublished  proprietary  data  and  has  no  effect  on  the  dynamic 
response  of  a  plate.  The  analytical  response  was  then  considered  valid  up  to 
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Che  point  at  which  a  failure  occurred  that  would  effect  the  dynamic 
response.  For  all  cases  considered,  this  failure  was  due  to  the  shear 
stress  in  the  top  ply  on  the  compressive  side  of  the  laminate.  In  this 
regard,  the  model  can  track  damage  to  a  limited  extent. 

It  must  be  noted  that  the  findings  of  these  analyses  are  for  a  specific 
impactor  and  may  not  hold  for  different  lmpactors  even  when  the  impact 
energy  is  the  same.  The  response  of  the  plate/impactor  system  depends  upon 
the  dynamics  of  both  parts.  Changing  the  mass  of  the  impactor  changes  its 
dynamics  and,  therefore,  changes  the  dynamics  of  the  system. 

To  demonstrate  this,  it  is  beneficial  to  consider  the  response  of  a 
system  in  generalized  coordinates.  With  this  formulation,  the  response  is 
given  by 


n 

u(x,y,z,t)  -  2  C  (t)  •  *.(x,y,z)  (43) 

i-1  1  1 

where  the  4>  are  eigenvectors  of  the  system,  and  the  C's  are  generalized 
coordinates.  The  first  eigenvector  characterizes  a  rigid  body  motion  of  the 
impactor,  the  second  is  the  first  mode  of  the  plate,  the  third  is  the  second 
mode  of  the  plate,  and  so  on. 

Time  histories  for  the  second  through  ninth  generalized  coordinates  of 
an  IM6/934  analysis  are  shown  in  Figure  20.  The  impact  scenario  for  this 
analysis  is  a  large  mass  traveling  at  a  low  velocity.  From  this  figure,  it 
can  be  seen  that  the  response  of  the  plate  is  primarily  a  first  mode 
response  with  a  small  second  mode  contribution.  In  contrast  to  this,  Figure 
21  shows  the  generalized  coordinates  for  the  same  plate  being  impacted  by  a 
small  mass  at  a  higher  velocity.  The  impact  events  for  these  analyses  are 
summarized  in  Table  6.  Note  that  the  impact  energy  is  the  same  for  both. 

There  are  many  differences  in  the  predicted  responses  for  these  two 
cases,  the  most  obvious  being  the  participation  of  higher  modes  than  the 
second.  The  initiation  of  damage  in  the  second  case  is  also  quite 
different.  As  in  the  analysis  for  the  large  impactor,  the  first  failure 
predicted  in  the  small  mass  analysis  is  matrix  cracking  in  the  tensile  side, 
followed  by  shear  failure  in  the  compressive  side.  In  the  first  analysis, 
the  matrix  cracking  is  predicted  at  about  500  microseconds  while  it  occurs 
around  25  microseconds  for  the  second.  Similarly,  the  first  analysis 
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predicted  the  shear  failure  at  about  13S0  microseconds  compared  to  35 
microseconds  for  the  second. 

More  important  than  the  times  at  which  these  failures  are  predicted, 
are  the  locations  at  which  they  occur.  The  large  mass  analysis  predicts  the 
shear  failure  in  the  upper  most  ply  of  the  laminate,  whereas  the  small  mass 
analysis  predicts  it  for  the  third  ply  in  from  the  top  surface.  Both 
analyses  predict  tensile  cracking  in  the  bottom  most  ply.  The  difference  in 
the  predicted  damage  location  means  that  the  spatial  stress  distribution  is 
dependent  upon  the  mass  of  the  impactor,  as  is  the  temporal  distribution. 

The  contribution  of  higher  modes  in  the  small  mass  analysis  is 
important  to  the  displacement  response  of  the  plate,  but  is  of  even  more 
importance  in  defining  the  stress  and  strain  response.  The  higher 
frequencies  are  associated  with  modes  characterized  by  shorter  wave  lengths 
resulting  in  higher  curvatures.  The  magnitude  of  the  generalized  coordinate 
multiplying  such  a  mode  may  then  be  relatively  small  and  still  have  a 
substantial  contribution  to  the  stress  distribution. 

For  example,  consider  the  bending  response  of  a  simply  supported 
uniform  beam.  Neglecting  shear  deformation  and  rotary  inertia,  the 
eigenvalues  for  this  beam  are  given  by 

Wj  -  (i  *2)  *  EI/mLA  (44) 

where  E  is  the  axial  modulus,  I  is  the  bending  moment  of  inertia,  m  is  the 
mass  of  the  beam,  and  L  is  its  length.  The  eigenvectors  associated  with 
these  eigenvalues  are 

w^(x)  -  2/mL  •  sin(i  w  x/L)  (45) 

From  (A3),  the  general  response  of  the  beam  is  then 


w(x,t)  +  Z  C.(t)  •  w  (x) 
i-1  1  1 


The  bending  stress  in  this  beam  is  given  by 


(A6) 
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(47) 


c  -  -Ez 
*  8x2 


Substituting  45  and  46  into  47  gives 


o  -  E  z  E  C.  (t)  {  2/mL  •  (i  */ L)2  •  sin(i  ir  x/L)  ]  (48) 

X  i-1  1 

Considering  the  contribution  of  the  first  two  symnetric  inodes  (i-1  and 
i-3)  to  the  total  stress  shows  the  importance  of  the  generalized  coordinate. 
The  contribution  from  the  first  mode  is  given  by 

o  -  E  z  C. (t)  2/mL  (*/L)2  sin(*  x/L)  (49) 

XI  1 


and  for  the  third  mode  by 


ax  -  9  E  z  C3(t)  2/mL  (x/L)2  sin(w  x/L)  (50) 


Comparing  (49)  to  (50)  it  can  be  seen  that  the  generalized  coordinate 
for  the  third  mode  can  be  nine  times  smaller  than  that  for  the  first  mode, 
and  still  have  a  contribution  to  the  stress  distribution  that  is  of  the  same 
order  of  magnitude  as  the  contribution  from  the  first.  The  comparison  with 
yet  higher  modes  becomes  even  more  striking  as  the  contribution  from  a  given 
mode  is  a  function  of  i  squared.  A  similar  expression  could  be  developed 
for  the  bending  modes  of  a  plate. 

The  generalized  coordinate  comparison  points  out  a  fundamental  problem 
in  understanding  low  velocity  impact.  A  large  number  of  variables  are 
associated  with  the  problem,  each  of  which  may  or  may  not  strongly  influence 
the  response  of  a  given  laminate.  It  is  doubtful  that  any  one  independent 
variable,  such  as  initial  impact  energy,  can  be  defined  that  will  be  a 
robust  measure  for  the  comparison  of  laminates.  Care  must  then  be  taken  to 
design  impact  tests  that  are  representative  of  the  expected  impact  scenario 
in  service. 
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It  nay  be  possible  that,  for  a  given  laainate,  an  inpact  test  nay  be 
devised  for  which  that  specific  laainate  perfoms  very  well.  A  nore 
informative  approach  to  the  problea  is  to  define  the  range  of  inpact 
scenarios  for  which  a  given  laainate  performs  well,  knowing  that  the 
definition  of  'performs  well'  is  also  somewhat  vague. 

The  analytical  model  developed  in  this  program  can  be  very  beneficial 
if  this  approach  is  used.  The  impactor  parameters,  most  notably  mass  and 
initial  velocity,  can  be  easily  varied,  allowing  the  analysis  of  a  wide 
range  of  impact  events.  These  analyses  can  then  be  used  to  identify  the  type 
of  application  for  which  this  laminate  is  suited.  As  with  any  analytical 
program,  a  number  of  select  cases  should  also  be  tested  experimentally  as  a 
check  on  the  numerical  results. 

A  number  of  extensions  to  the  existing  model  are  possible  and  should  be 
considered  for  further  development.  Most  important  is  the  use  of  a  non¬ 
linear  formulation  for  the  analysis  of  thin  plates  followed  by  including 
mechanisms  for  tracking  damage.  As  previously  mentioned,  this  first 
extension  is  well  defined  and  poses  no  Insurmountable  analytical  problem. 

It  does,  however,  require  the  use  of  the  fastest  computer  available. 

The  second  extension  is  not  so  clear  cut,  and  does  pose  some 
significant  analytical  challenges.  Modifying  ply  properties  to  account  for 
matrix  cracking  or  even  fibers  breaking  is  fairly  routine,  but  predicting 
the  growth  of  a  crack  is  not.  Finite  element  methods  for  predicting  the 
growth  of  a  crack  have  been  developed  (Reference  15  and  16)  ,  but  are  quite 
complex  even  for  static  loads.  The  nature  of  the  cracks  that  develop  due  to 
low  velocity  impact  further  complicate  the  issue  since  they  are  trans-ply 
rather  than  inter-ply  cracks. 

The  accumulation  of  damage  in  some  laminates  may  occur  in  almost 
discrete  steps  with  a  damage  zone  being  initiated  and  fully  developed  over  a 
very  short  period  of  time.  An  approximate  model  which  exploits  this  effect 
may  be  postulated  for  tracking  damage  in  these  laminates.  In  this  model  the 
initial  stiffness  of  the  laminate  would  be  used  to  predict  the  response  up 
to  the  initiation  of  a  shear  crack.  A  fracture  mechanics  model  would  then 
be  used  to  predict  the  growth  of  the  crack  up  to  the  point  at  which  it 
became  stable.  A  new  laminate  stiffness  would  be  formulated  which  included 
the  damage  zone,  and  used  to  predict  the  dynamic  response  up  to  the  next 
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damage  initiation.  Obviously,  this  is  a  simplified  approach  to  a  very 
complex  problem  but  may  warrant  a  feasibility  study. 

It  is  recommended  that  any  further  extensions  of  the  existing  model  be 
made  on  the  two-dimensional  case.  It  is  not  obvious  that  results  predicted 
by  the  three-dimensional  model  will  yield  any  information  about  a  laminate 
that  is  of  a  more  general  nature  than  2-D  results.  Of  course  there  are 
problems  associated  with  the  testing  of  beams,  such  as  edge  effects,  but  the 
reduced  analytical  effort  is  very  attractive.  At  the  very  least,  the  2-D 
model  could  be  extended  for  performing  feasibility  studies  before  expending 
the  resources  on  the  complete  3-D  formulation. 
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Table  1.  AS4/3502  Analysis  Summary 


Analysis 

Plies* 

Support 

Diameter 

(mm) 

Impactor 

Mass 

(kg) 

Iapactor 

Velocity 

(m/s) 

1 

16 

50 

2.49 

1.56 

2 

16 

75 

5.04 

1.56 

3 

8 

25 

2.49 

1  56 

4 

8 

50 

2.49 

1.56 

5 

8 

75 

2.49 

1.56 

*  8  ply  laminate  -  [0 ,+45 , -45 , 90] s 

16  ply  laminate  -  { 0 , +45 , -45 , 90] 2s 
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Table  2.  Material  Properties  for  AS4/3502  Analyses 


Unidirectional  Ply  Properties: 


E1 

-  130.3  GPa 

E2 

-  9.65  GPa 

*12 

0.35 

*23 

0.40 

G12 

5.93  GPa 

P 

-  4.08xl0'3 

t 

-  0 . 13  mm 

sT 

*1 

-  1793.  MPa 

sG 

-  1724.  MPa 

X 

A 

-  64 .  MPa 

A 

-  138.  MPa 

S12 

-  109 .  MPa 

m 

CM 

w 

69 .  MPa 

Neat 

Resin  Propertie 

3.79  GPa 
0.36 

3.27x10  3g/cc 
38.  MPa 

69.  MPa 

62.  MPa 
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Table  3. 


Contact  Law  Parameters  for  AS4/3502  Analyses 


K  -  1.07xl05N/mm3/2 

n  -  1.5 

m  -  2.5 

Acr  -  0.08  mm 

1  -  0.4 
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Analysis 


1 


Table  4.  IM6/934  Analysis  Summary 


Plies  Support  Impactor 

(mm)  Mass 

(kg) 

32  76.2  X  127.0  3.17 


*  [0,90,+45,-45]4s 


Impactor 

Velocity 

(m/s) 

1.10 
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Table  5.  Aspect  Ratio  Summary 


Laminate 

Thickness 

(mm) 

Support 

(mm) 

Aspect  Ratio 

AS4/3502 

2.1133  (16  ply) 

50  (Dia.) 

25 

AS4/3502 

2.1133  (16  ply) 

75  (Dia.) 

35 

AS4/3502 

1.0566  (8  ply) 

25  (Dia.) 

24 

AS4/3502 

1.0566  (8  ply) 

50  (Dia.) 

47 

AS4/3502 

1.0566  (8  ply) 

75  (Dia.) 

71 

1M6/934 

4.1453  (32  ply) 

76.2  X  127.0 

18  X  31 
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Table  6.  Generalized  Coordinate  Analysis  Summary  (IM6/934) 


Analysis 

Support 

Impactor 

Impactor 

Impact 

(mm) 

Mass 

Velocity 

Energy 

(kg) 

(m/s) 

(N-m) 

1 

76.2  X  127.0 

3.17 

1.56 

3.84 

2 

76.2  X  127.0 

0.016 

21.77 

3.84 
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Figure  1.  Experimental  Contact  Force  History  Comparison  for  8 
Ply  Laminates  With  75  mm  Diameter  Support 
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Experimental  Contact  Force  History  Comparison  for  16  Piy 
Laminates  With  75  mm  Diameter  Support 


* 


1.  Model  Concept 


3.  Reduced  Substructures 
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Figure  4.  2-D  Element  Nomenclature 
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Figure  5.  Special  3-D  Geometry 
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Figure  6.  General  3-D  Geometry 
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Figure  8.  Typical  Force-Deflection  Diagram  for  the  Contact  Spring 
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Figure  13.  Normalized  Stress  Histories  for  16  Ply  AS4/3502 
Laminate  with  50  mm  Diameter  Support 
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Figure  14.  Contact  Force  History  for  IM6/934  Laminate.  Impact 
Energy  =  1.92  N-m. 
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19.  Extensional  Strain  Comparison  for  8  Ply  AS4/3502 
Laminate  with  50  mm  Diameter  Support 
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APPENDIX  -  SOFTWARE  DESCRIPTION 


OVERVIEW 

The  three  programs  LAMCEL,  1AMSUB  and  TRADYR  constitute  an  analysis 
system  for  predicting  the  transient  dynamic  response  of  laminated 
composites,  and  in  particular,  the  initiation  of  damage  due  to  low  velocity 
impact.  The  dynamic  analysis  utilizes  the  finite  element  method  and 
consists  of  three  stages;  mass  and  stiffness  characterization,  numerical 
integration  of  the  equations  of  motion,  and  postprocessing  for  stresses  and 
strains.  Either  LAMCEL  or  LAMSUB  may  be  used  for  the  first  and  third 
stages,  and  TRADYR  is  used  for  the  second. 

LAMCEL  is  used  to  model  laminated  beams  in  two - dimens ions ,  or 
rectangular  plates  in  three.  This  program  exploits  the  regular  geometry  of 
the  uniform  beam  or  rectangular  plate  to  optimize  the  generation  of  the  mass 
and  stiffness  of  the  structure.  The  dynamic  properties  of  a  representative 
'cell'  are  generated  by  LAMCEL,  and  used  to  define  the  properties  of  the 
entire  beam  or  plate  in  TRADYR. 

LAMSUB  may  only  be  used  to  model  three-dimensional  plates,  but  is  not 
restricted  to  rectangular  geometries.  This  program  generates  the  mass  and 
stiffness  matrices  for  any  arbitrary  sub-region  of  the  plate. 

TRADYR  combines  the  mass  and  stiffness  matrices  generated  by  LAMCEL  or 
LAMSUB  to  produce  the  mass  and  stiffness  matrix  for  the  beam  or  plate  to  be 
analyzed.  To  this  it  adds  the  dynamic  properties  of  the  impactor  resulting 
in  the  equations  of  motion  for  the  complete  plate/impactor  system.  These 
equations  of  motion  are  then  integrated  numerically  for  the  displacement 
history  of  the  system. 

The  displacement  history  generated  by  TRADYR  is  used  by  either  LAMCEL 
or  LAMSUB  to  calculate  the  stress  and  strain  histories  of  the  structure.  A 
failure  analysis  is  also  performed  by  these  post  processors  to  determine  the 
onset  of  damage.  It  is  assumed  that  the  plate  behaves  as  a  linear  elastic 
structure  up  to  the  initiation  of  damage  with  all  nonlinearities  being 
included  in  the  contact  model. 
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ANALYSIS  PROCEDURE 


Extensive  use  of  substructuring  is  used  in  the  analysis,  resulting  in 
models  that  include  adequate  detail  without  being  prohibitively  large.  With 
this  technique  a  small  portion  of  the  structure  (substructure)  is  modeled  in 
great  detail  with  a  large  number  of  degrees  of  freedom.  The  equations  of 
motion  for  this  substructure  are  then  reduced  to  include  a  much  smaller 
number  of  'master'  degrees  of  freedom.  All  other  degrees  of  freedom 
(slaves)  will  be  expressed  in  terms  of  the  masters. 

Several  substructures  may  then  be  combined  to  provide  a  coarse  model 
for  the  dynamic  analysis.  By  retaining  all  degrees  of  freedom  along  the 
boundaries  between  substructures  as  masters,  a  model  having  displacement 
compatability  will  be  preserved.  A  further  reduction  In  the  assembled 
dynamic  model  is  achieved  by  selecting  a  number  of  'governor'  degrees  of 
freedom  from  the  masters.  The  model  is  then  condensed  by  defining  the 
remaining  master  degrees  of  freedom  in  terms  of  the  governors. 

In  a  uniform  beam  or  rectangular  plate  model,  the  structure  may  be 
divided  into  a  number  of  uniform  rectangular  cells.  These  cells  are 
identical,  allowing  the  stiffness  and  inertial  properties  of  a  single  cell 
to  be  generated  once  and  used  for  all.  For  plates  of  arbitrary  geometry, 
definition  of  a  single  representative  cell  may  not  be  possible.  These 
plates  must  be  modeled  using  a  number  of  unique  substructures. 

To  complete  the  impact  model  a  lumped  mass  representing  the  impactor  is 
tied  to  a  governor  degree  of  freedom  on  the  plate  throught  a  non¬ 
conservative,  non-linear  spring.  The  plate  then  sees  the  impactor  as  simply 
a  time  varying  force.  The  magnitude  of  the  force  and  the  amount  of  energy 
dissipated  being  functions  of  the  relative  motion  between  the  plate  and 
impactor . 

The  equations  of  motion  for  the  resulting  dynamic  system  are  integrated 
to  yield  displacements,  velocities  and  accelerations  as  functions  of  time. 
The  integration  is  performed  numerically,  and  displacement  histories  are 
stored  on  disk  for  subsequent  post-processing.  A  restart  capability  is 
included  in  the  numerical  integration  routine,  allowing  the  analysis  to  be 
stoped,  results  reviewed,  and  the  analysis  restarted. 

Displacement  solutions  predicted  by  the  reduced  dynamic  model  may  be 
used  to  define  the  boundary  conditions  for  refined  stress  analyses  of  the 
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original  substructures.  The  substructure  concept  allows  the  stress  analysis 
to  be  performed  only  for  the  regions  of  interest,  such  as  in  the  vicinity  of 
the  impactor.  To  further  minimize  the  computational  effort,  elements  within 
a  given  substructure  are  specified  for  subsequent  post  processing.  Stresses 
and  strains  are  only  calculated  for  those  elements  which  Lave  been 
specified. 

Quantities  which  may  be  calculated  for  each  element  are  stresses  and 
strains  in  global  and/or  material  coordinates.  Global  strains  are 
calculated  first  and  transformed  into  material  coordinates  for  subsequent 
calculations.  Stresses  and  strains  are  volume  averaged  over  each  ply  within 
an  element,  and  area  averaged  over  interfaces.  Unique  interface  properties 
may  be  defined  for  stress  calculation  purposes.  The  interface  is  assumed  to 
have  no  volume  and,  therefore,  does  not  contribute  to  the  stiffness  of  the 
element . 

A  failure  analysis  is  performed  for  each  ply  in  each  selected  element, 
and  for  each  interface.  The  model  does  not  modify  the  stiffness  of  the 
structure,  based  on  the  prediction  of  failure,  and  may  not  be  valid  for 
analysis  past  the  point  of  initial  damage. 

ELEMENT  DESCRIPTION 

Two-  or  three-dimensional  solid  finite  elements  are  used  to  define  the 
mass  and  stiffness  matrices  of  the  beam  or  plate  respectively.  The  geometry 
of  the  two-dimensional  element  is  defined  by  four  comer  nodes  (Figure  Al), 
and  the  three-dimensional  element  by  eight  comer  nodes  (Figures  A2  and  A3). 
Elements  are  allowed  to  have  more  than  one  ply  throught  the  thickness  but 
must  contain  complete  plies. 

Nodal  variables  for  these  elements  consist  of  translations  in  each 
coordinate  direction  and  the  partial  derivative  of  each  translation  with 
respect  to  each  coordinate.  This  results  in  six  degrees  of  freedom  per  node 
for  the  two-dimensional  element,  and  twelve  for  the  three-dimensional 
element.  Cubic  polynomials  are  used  as  shape  functions  to  define  the 
translational  displacement  fields  for  these  elements. 

Throughout  this  program,  nodal  degrees  of  freedom  are  referred  to  by 
the  nomenclature  Unm,  where  n  is  an  integer  from  1  to  3,  and  m  is  an  integer 
from  0  to  3.  Unm  is  then  interpreted  as  the  partial  derivative  of  the 
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Figure  A-l.  2-D  Element  Description 


Eight  Element  Cell 


Node  (NLAYER+1) *n 
Node  (NLAYER+1) *n-l 


Node  Column 

<  i 

<>  Node  (NLAYER+1) *n-NLAYER 
Figure  A-3.  3-D  General  Element  Description 


translation  in  the  n  direction,  with  respect  to  the  a  coordinate.  The  only 
exception  is  when  a  is  0.  For  this  case,  UnO  is  interpreted  as  the 
translation  in  the  n  direction.  Relative  to  displacenents,  direction  1  is 
the  global  X  direction,  2  is  Y  and  3  is  Z,  where  the  X,  Y  and  Z  axes  form  a 
right  handed  coordinate  system.  The  Z  direction  is  always  assumed  to  be 
through  the  thickness  of  the  laminate,  and  two-dimensional  models  are  always 
defined  in  the  X-Z  plane. 

Consistent  mass  matrices  are  generated  using  linear  'Serindipity'  shape 
functions.  The  resulting  matrix  only  associates  mass  with  translational 
degrees  of  freedom.  The  algorithm  used  to  numerically  integrate  the 
equations  of  motion  requires  that  the  reduced  mass  matrix  be  invertable. 
Since  there  is  no  mass  associated  with  degrees  of  freedom  other  than 
translations,  the  reduced  mass  matrix  for  a  system  which  included  partial 
derivatives  as  governor  degrees  of  freedom  would  have  zeros  on  the  diagonal. 
Such  a  matrix  would  not  be  invertable,  disallowing  the  inclusion  of  partial 
derivatives  as  governor  degrees  of  freedom. 

LAMCEL  AND  LAMSUB  -  PROGRAM  DESCRIPTIONS  AND  INPUT  FORMAT 

The  programs  LAMCEL  and  LAMSUB  are  similar  in  that  both  are  used  to 
generate  the  mass  and  stiffness  matrices  for  a  sub-region  of  a  plate,  and  to 
perform  stress  analyses  for  those  same  sub-regions.  Except  for  the 
definition  of  the  finite  element  mesh  used  by  each  program,  and  minor 
differences  in  the  control  parameters,  data  input  to  the  two  programs  is 
identical . 

LAMCEL  is  used  to  generate  the  mass  and  stiffness  matrices  for  a  three- 
dimensional  cell  that  is  representative  of  the  mass  and  stiffness  of  a  small 
region  of  a  rectangular  plate.  In  addition,  LAMCEL  may  be  used  to  define  a 
two-dimensional  cell  for  analyzing  the  impact  response  of  a  uniform  beam. 
Two-dimensional  cells  may  assume  states  of  either  plane  stress,  or  plane 
strain. 

LAMSUB  can  only  be  used  to  analyze  plates  in  three-dimensions  but 
places  no  restriction  on  the  geometry  of  the  plate.  The  penalty  paid  for 
this  increased  flexiblility  is  increased  run  time,  and  more  involved  model 
definition. 
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For  generation  analyses,  both  programs  read  all  input  from  an  ASCII 
file  named  FILE01.DAT.  Output  is  written  to  an  ASCII  file  named  FILE02.DAT, 
and  a  number  of  binary  files  named  by  the  user.  These  binary  files  are  used 
in  subsequent  dynamic  analyses  by  TRADYR,  and  back  substitution  analyses  by 
either  LAMCEL  or  LAMSUB,  as  shown  in  Figures  A-4  and  A* 5. 

Back  substitution  analyses  also  require  input  from  FILE01.DAT,  and 
write  to  FILE02.DAT.  In  addition,  a  binary  input  file  containing  the 
displacement  history  predicted  by  TRADYR  must  also  be  furnished.  The  binary 
files  written  during  the  cell  generation  pass  are  also  utilized  during  back 
substitution. 

Data  is  read  sequentially  from  FILE01.DAT  with  commas  separating  data 
items  within  a  line  of  input.  Successive  commas  will  be  interpreted  as  0  or 
0.0  for  integer  or  real  input.  Detailed  descriptions  of  the  input  data  are 
included  in  the  following  sections. 

A.  Control 

TITLE  78  character  alphanumeric  title. 

NCHECK , NATYPE , NORDER , IMPRNT  The  data  check  flag,  analysis 

type  flag,  numerical  integration 
order,  and  the  expanded  material 
property  printout  flag. 

NCHECK  -  0  (default),  data  check  only. 

-  1,  analysis 

In  LAMCEL,  NATYPE  -  0,  back  substitution. 

-  1,  plane  stress  generation. 

-  2,  plane  strain  generation. 

-  3,  3-D  generation. 

In  LAMSUB,  NATYPE  -  0,  back  substitution. 

-  1,  3-D  generation. 

2  <  NORDER  <  4 

IMPRNT  -  0  (default),  print  matrices. 

-  1,  do  not  print  matrices. 
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Generation  Solution  Back  Substitution 


Fiqure  A-4.  T.AMCFL  File  Ueaqe 


Solution  Back  Substitution 


Figure  A-5.  LAMSUB  File  Usage 


NAME1 


Back  substitution  data  storage 
file  name.  (40  character  maximum) 


NAME2  Lamina  definition  storage  file 

name.  (40  character  maximum) 

NAME3  Element  definition  storage  file 

name.  (40  character  maximum) 

NAME4  Constraint  definition  storage 

file  name.  (40  character  maximum) 

NAME5  Strain  recovery  data  storage  file 

name.  (40  character  maximum) 

NAME6  If  NATYPE  >  0,  constrained, 

condensed  cell  mass  and  stiffness 
matrix  storage  file  name.  If 
NATYPE  -  0,  displacement  vector 
storage  file  name.  (40  character 
maximum) 

The  program  will  branch  at  this  point  depending  upon  the  value  of 
NATYPE.  If  NATYPE  is  set  to  zero,  the  program  will  expect  back  substitution 
(Section  B)  information  next.  If  NATYPE  is  set  to  any  other  valid  value  it 
will  branch  to  the  laminate  description  section  (Section  C) . 

NORDER  specifies  the  order  of  the  numerical  integration  algorithm  used 
to  calculate  the  stiffness  matrices  of  elements.  Two,  three  or  four  point 
integration  may  be  specified.  Expanded  material  property  printout  consists 
of  the  stress-strain  matrix  for  each  material  in  material  coordinates,  and 
the  stress-strain  matrix  for  each  lamina  in  global  coordinates.  For  back 
substitution  analyses  (NATYPE-0) ,  NORDER  and  IKPRNT  are  not  needed. 
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B.  Back  Substitution  Information 


NDISP.NSUBS 


IDPRNT,(NPRINT< I) ,1- 


NEBS 


IELEM 


The  index  number  of  the 
displacement  set  to  be  used,  and 
the  index  number  of  the 
substructure  (cell)  to  be 
analyzed. 

,7)  The  expanded  displacement  vector 

printout  flag  and  the  stress  and 
strain  printout  flags.  When  a 
flag  is  set  to  0  (default) ,  the 
associated  data  set  will  be 
printed,  if  set  to  1,  the  data 
will  not  be  printed. 

IDPRNT  -  displacement  vector. 

NPRINT(l)  -  global  stresses. 

NPRINT(2)  -  global  strains. 

NPRINT(3)  -  local  stresses. 

NPRINT(4)  -  local  strains. 

NPRINT(5)  -  ply  printout. 

NPRINT(6)  -  interface  printout. 

NPRINT(7)  -  surface  printout. 

The  number  of  elements  for  which 
stresses  and  strains  will  be 
calculated.  The  program  will  look 
for  NEBS  element  numbers  following 
this  entry. 

The  index  number  of  the  element 
for  which  stresses  and  strains 
will  be  calculated. 
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After  performing  a  stress  analysis  on  the  last  element  specified  fh 
this  set,  the  program  will  loop  back  to  input  the  next  values  for  HDISP  and 
NSUBS .  This  loop  will  continue  until  an  end  of  file  is  encountered  while 
trying  to  read  NDISP. 

The  files  specified  by  NAME1  through  NAME5  in  the  control  section  for  a 
back  substitution  analysis  must  be  the  same  as  those  defined  in  the 
generation  analysis.  The  file  specified  by  NAME6  is  the  same  as  that 
defined  by  NAME1  in  TRADYR  (see  Figures  A-4  and  A-5).  This  file  contains 
the  displacement  vector  for  each  substructure  in  the  assembled  model  at  each 
print  step.  NDISP  is  the  number  of  the  print  step  for  which  the  back 
substitution  is  being  performed,  and  NSUBS  is  the  index  number  of  the 
substructure  being  analyzed. 

Stresses  and/or  strains  may  be  calculated  in  global  and/or  local 
coordinates.  The  quantities  may  be  calculated  at  any  or  all  of  three 
locations;  volume  averages  within  plies,  area  averages  of  ply  interfaces, 
and  area  averages  at  the  surfaces  of  the  laminate.  The  quantities 
calculated,  and  the  locations  at  which  they  are  calculated,  are  controlled 
by  the  flags  NPRINT(2)  through  NPRINT(7). 


C.  Laminate  Description 

NMAT 


LABEL 


El ,E2 ,NU12 , NU23.G12 ,RHO 


S1T,S1C,S2T,S2C,S12,S23 


The  number  of  materials  used  in 
the  laminate.  The  program  will 
look  for  NMAT  sets  of  material 
property  input  following  this 
entry . 

40  character  alphanumeric 
material  label. 

Elastic  constants  for  this 
material . 

Strengths  for  this  material. 
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NPLY , IFMAT 


The  number  of  plies  in  Che 
laminate,  and  the  index  number  of 
the  material  to  be  used  for 
calculating  interface  stresses 
and  strains.  The  program  will 
look  for  NPLY  lamina  definitions 
following  this  entry. 

MAT, T, THETA  The  material  index  number, 

thickness  and  orientation  of  the 
ply. 

Material  property  sets  are  assigned  sequential  index  numbers  based  on 
the  order  in  which  they  are  input.  Each  set  defines  the  properties  for  a 
unidirectional  composite  material  for  subsequent  use  in  defining  a 
particular  laminate.  Elastic  constants  El  and  E2  are  the  axial  and 
transverse  elastic  moduli,  NU12  and  NU23  are  the  major  Possion's  ratios  in 
the  1-2  and  2-3  planes,  G12  is  the  shear  mudulus  in  the  1-2  plane,  and  RHO 
is  the  density  in  mass  density  units.  All  other  elastic  constants  are 
calculated  assuming  that  the  material  is  transversely  isotropic. 

Strengths  SIT  and  SIC  are  the  tensile  and  compressive  axial  strengths, 
S2T  and  S2C  are  the  tensile  and  compressive  transverse  strengths,  and  S12 
and  S23  are  the  shear  strengths  in  the  1-2  and  2-3  planes.  As  with  elastic 
constants  all  other  strengths  are  assigned  assuming  that  the  material  is 
transversely  isotropic  in  strength. 

The  laminate  is  defined  by  specifying  the  type  of  material  and 
orientation  of  each  ply  starting  at  the  bottom  of  the  laminate  and  working 
up.  Index  numbers  are  assigned  to  the  lamina  sequentially  based  on  the 
order  in  which  they  are  input.  The  index  numbers  are  subsequently  used  to 
specify  which  lamina  are  contained  in  a  given  element.  The  lamina 
definition  parameter  MAT  specifies  which  set  of  material  properties  to  use 
for  a  given  lamina  by  referencing  that  materials  index  number.  The 
parameter  THETA  specifies  the  rotation  of  the  material  1  axis  in  the  global 
X-Y  plane.  THETA  is  measured  from  the  global  X  axis  to  the  material  1  axis, 
in  degrees,  positive  by  the  right  hand  rule. 
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All  interfaces  are  assumed  to  have  the  same  elastic  properties  and 
strengths.  The  properties  specified  by  the  parameter  IFMAT  are  only  used 
for  stress  and  strain  calculations,  and  do  not  contribute  to  the  stiffness 
of  the  laminate. 

Following  the  laminate  description  input,  the  program  will  look  for 
cell  (Section  D.l)  or  substructure  (Section  D.2)  definition  input,  depending 
upon  which  program  is  being  run. 

D.l  LAMCEL  Cell  Definition 

NE  The  number  of  elements  in  the 

cell.  The  program  will  look  for 
NE  sets  of  element  definition 
input  following  this  entry. 

ELEN, EWID, IPLYB, IPLYT  The  length  and  width  of  the 

element,  and  the  index  numbers  of 
the  bottom  and  top  plies  in  the 
element . 

(NODE(I) , 1-1 ,NNE)  The  index  numbers  for  the  comers 

of  the  element. 

If  NATYPE  -  1  or  2,  NNE  -  4 
3,  NNE  -  8 

A  cell  is  a  refined  finite  element  model  of  a  small  rectangular  region 
of  a  laminate.  Each  element  in  the  cell  is  assumed  to  be  rectangular,  and 
is  defined  by  a  length  (ELEN)  and  width  (EWID)  in  the  X-Y  plane,  and  the 
plies  through  its  thickness  (IPLYB  and  IPLYT).  The  thickness  of  an  element 
is  defined  by  the  sum  of  the  thicknesses  of  the  plies  contained  in  it,  and 
is  never  specified  explicitly.  For  two-dimensional  analyses  (NATYPE  -  1  or 
2),  EWID  may  be  omitted,  and  will  default  to  1. 

Node  positions  are  not  specified  explicitly,  but  are  defined  by  the 
element  geometries  and  connectivity.  Element  node  numbers  must  be  input  in 
the  order  show  in  Figure  A-l  for  two-dimensional  analyses,  and  Figure  A- 2 
for  three-dimensional  analyses. 
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Elements  are  assigned  sequential  index  numbers  based  on  the  order  in 
which  they  are  input.  These  index  numbers  serve  as  a  means  for  referencing 
a  particular  element  in  subsequent  post  processing.  The  program  will  expect 
constraining  information  (Section  E)  following  the  cell  definition  input. 

D.2  LAMSUB  Substructure  Definition 


NLAYER 


1PLYB.IPLYT 


NC 


NCOL.X.Y.BETA 


NS 


NCOL(I) ,1-1,4 


The  number  of  layers  ot  elements 
through  the  thickness  of  the 
laminate.  The  program  will  look 
for  NLAYER  layer  definitions 
following  this  entry. 

Index  numbers  for  the  bottom  and 
top  plies  in  the  layer. 

The  number  of  columns  of  nodes  to 
be  defined.  The  program  will 
look  for  NC  column  definitions 
following  this  entry. 

The  index  number  of  the  column, 
its  X  and  Y  position,  and  the 
rotation  of  the  nodal  coordinate 
system  relative  to  the  global 
axes.  Beta  is  input  in  degrees. 

The  number  of  element  stacks  to 
be  defined.  The  program  will 
look  for  NS  stack  definitions 
following  this  entry. 

The  node  columns  defining  the 
comers  of  this  stack. 
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The  finite  element  aesh  for  the  three-dimensional  substructure  is 
defined  by  a  two-dimensional  mesh  of  quadrilaterals  in  the  X-Y  plane,  and 
information  about  the  number  and  thickness  of  elements  through  the  thickness 
of  the  laminate.  Each  of  the  quadrilaterals  represents  a  stack  of  NLAYER 
elements  having  uniform  geometry  in  the  Z  direction.  Hie  corners  of  the 
quadrilateral  represent  the  position  of  a  column  of  nodes  which  all  have  the 
same  X  and  Y  coordinates  but  different  Z  positions. 

Each  layer  consists  of  one  or  more  unidirectional  plies  and  represents 
a  one  element  thick  mesh  in  the  X-Y  plane.  The  thickness  of  a  layer  is 
defined  by  the  sum  of  the  thicknesses  of  the  plies  contained  in  it,  and  is 
never  specified  explicitly.  Layer  interfaces  and  ply  interfaces  must  not  be 
confused.  Layer  interfaces  must  occur  at  ply  interfaces,  but  ply  interfaces 
do  not  necessarily  have  to  be  layer  interfaces.  The  bottom  most  ply  in  the 
first  layer  must  be  ply  1  (IPLYB(l)-l) ,  and  subsequent  layers  must  be  input 
in  order  of  increasing  ply  index.  The  index  number  of  the  bottom  ply  of 
layer  n  ( IPLYB(n) )  must  be  one  greater  than  the  index  number  of  the  top  ply 
of  layer  n-1  ( I PLYT ( n - 1 ) ) . 

The  program  automatically  places  nodes  through  the  thickness  of  the 
laminate  for  each  location  at  which  a  nodal  column  in  defined.  The  first  of 
these  nodes  is  placed  at  the  bottom  of  the  first  layer;  nodes  are  then 
located  at  each  layer  interface,  starting  with  the  bottom  layer  and  working 
toward  the  top;  the  last  node  in  the  column  is  placed  at  the  top  of  the  last 
layer.  Each  column  then  consists  of  NLAYER+1  nodes.  Each  node  is  assigned  a 
unique  index  number,  based  on  the  index  number  of  the  column  it  is  located 
in.  The  index  number  of  the  first  node  in  column  n  is  given  by 
(KLAYER+l)*(n-l)+l .  Index  numbers  for  subsequent  nodes  in  that  column  are 
then  incremented  by  1. 

The  shape  of  an  element  stack  may  be  any  arbitrary  quadrilateral,  in 
the  X-Y  plane,  and  is  uniform  in  the  Z  direction.  The  edges  of  the  stack 
are  defined  by  the  index  numbers  of  the  four  node  columns  at  the  comers  of 
the  quadrilateral.  The  order  in  which  these  indices  are  input  must  follow 
the  convention  shown  in  Figure  A- 3.  Each  stack  consists  of  NLAYER  elements 
through  the  thi -kness  of  the  laminate.  Elements  are  assigned  sequential 
index  numbers  based  on  the  order  in  which  the  stacks  are  input.  Within  a 
stack  the  element  numbering  begins  with  the  bottom  layer,  incrementing  by  1 
as  the  top  layer  Is  approached.  The  index  number  for  the  nth  element 


A- 17 


(layer)  In  the  mth  stack  of  a  laminate  with  NLAYER  layers  Is  given  by 
NLAYEP*(m-l)+n.  These  Index  numbers  serve  as  a  means  for  referencing  a 
particular  element  in  subsequent  post  processing. 

The  program  will  expect  constraint  Information  (Section  E)  following 
the  cell  definition  Input. 


E.  Constraints 


NZDF 


NODE .DIR 


NCS 


NTS , DIR 


NODE(I) ,1-l.NIS 


The  number  of  constrained 
(fixed)  displacements.  The 
program  will  look  for  NZDF 
constraint  definitions  following 
this  entry. 

The  node  which  is  constrained, 
and  the  direction  in  which  it  is 
constrained. 

The  number  of  coupled  sets,  the 
program  will  look  for  NCS  sets 
of  coupled  set  input  following 
this  entry 

The  number  of  nodes  in  the  set, 
and  the  direction  in  which  they 
are  coupled. 

Index  numbers  of  the  coupled 
nodes,  input  in  groups  of  10  per 
line . 


NCE  The  number  of  constraint 

equations.  The  program  will 
look  for  NCE  sets  of  constraint 
equation  input  following  this 
entry. 
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NTE  The  number  of  terns  in  the 

equation.  The  program  will  look 
for  NTE  constraint  equation 
terms  following  this  entry. 

NODE, DIR, VALUE  Constraint  equation  term. 

Three  types  of  constraints  may  be  defined  for  a  given  degree  of 
freedom.  The  degree  of  freedom  may  be  specified  as  having  zero  magnitude, 
the  same  magnitude  as  a  similar  degree  of  freedom,  or  a  magnitude  defined  by 
a  linear  combination  of  several  other  degrees  of  freedom.  These  types  of 
constraints  are  referred  to  as  fixed  displacements ,  coupled  sets,  and 
constraint  equations,  respectively. 

The  constraint  definition  NODE, DIR  specifies  that  there  will  be  no 
motion  of  at  node  NODE  in  the  DIR  direction.  This  degree  of  freedom  is  then 
removed  from  the  problem  and  may  not  be  used  in  any  subseqent  constraint 
definition  or  master  degree  of  freedom  specifications.  If  NODE  is  input  as 
zero,  all  nodes  will  be  constrained  in  the  DIR  direction. 

Coupled  sets  force  all  of  the  specified  nodes  to  have  the  same  motion 
in  the  direction  specified  by  DIR.  The  first  node  specified  in  the  set 
remains  active  in  that  direction.  For  all  other  nodes  In  the  set,  that 
degree  of  freedom  will  be  constrained  out  of  the  problem  and  may  not  be  used 
in  subsequent  constraint  definitions  or  master  degree  of  freedom 
specifications . 

Constraint  equations  are  linear  equations  which  define  the  motion  of 
one  degree  of  freedom  in  terms  of  one  or  more  other  degrees  of  freedom.  The 
terms  in  the  constraint  equation  input  set  are  interpreted  as 

VALUE (1)*[N0DE(1) ,DIR(1) J  -  VAUJE(2)*[NODE(2) ,DIR(2) ] 

+  VALUE(3)*[NODE(3) ,DIR(3) ] 


+  VALUE ( NTE) *[NODE(NTE) ,DIR(NTE) ] 
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vbere  [NODE(I) ,DIR(I) ]  specify  the  degree  of  freedon  associated  with  the  Ith 
tens,  and  VALUE(I)  is  the  coefficient  multiplying  that  term.  The  degree  of 
freedom  specified  in  the  first  term  of  the  set  is  constrained  out  of  the 
problem,  and  may  not  be  used  in  subsequent  constraint  definitions  or  master 
degree  of  freedom  specifications. 

Direction  specifications  (DIR)  are  three  character  labels  of  the  form 
Unm,  where  n  is  an  integer  from  1  to  3 ,  and  m  is  an  integer  from  0  to  3. 

When  m  is  zero,  UnO  is  interpreted  as  the  translation  in  the  n  direction. 
When  m  is  not  zero,  Unm  is  interpreted  as  the  partial  derivative  of  the 
translation  in  the  n  direction  with  respect  to  the  m  coordinate.  For  LAMCEL 
analyses,  direction  1  is  the  global  X  direction,  2  is  Y,  and  3  is  Z.  For 
LAMSUB  analyses,  the  nodal  coordinate  system  (1,2,3)  may  be  rotated  relative 
to  the  global  axes  (X,Y,Z).  Constraints  are  always  specified  in  the  nodal 
coordinate  system.  Directions  1  and  3  are  the  only  valid  directions  for 
two-dimensional  LAMCEL  analyses. 

When  a  model  is  generated  that  does  not  utilize  a  particular  type  of 
constraint  the  parameter  which  defines  the  number  of  sets  of  data  for  that 
type  of  constraint  (NZDF,  NCS  or  NCE)  must  still  be  included  as  zero. 
Following  the  constraint  input,  the  program  will  look  for  master  degree  of 
freedom  specifications  (Section  F) . 

F.  Master  Degrees  of  Freedom 

NMDF  The  number  of  master  degrees  of 

freedom.  The  program  will  look 
for  NMDF  master  degree  of 
freedom  specifications  following 
this  entry. 

NODE, DIR  The  node  number  and  direction 

that  define  the  master  degree  of 
freedom. 

Master  degrees  of  freedom  are  assigned  index  numbers  sequentially  based 
on  the  order  in  which  they  are  Input.  The  first  definition  input  will 
become  master  d.o.f.  1  of  the  cell  regardless  of  the  values  of  NODE  and  DIR; 
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the  second  will  become  master  d.o.f.  2,  and  so  on.  Rows  and  columns  of  the 
constrained,  condensed  mass,  and  stiffness  matrices  for  the  cell  will  be 
arranged  to  reflect  this  numbering  system. 

The  direction  parameter  DIR  may  be  any  valid  label  as  defined  in  the 
constraint  section  (Section  E) .  In  addition,  DIR  may  be  assigned  the  label 
ATI.,  in  which  case  all  unconstrained  degrees  of  freedom  at  the  specified 
node  will  we  retained  as  master  degrees  of  freedom.  If  NODE  is  input  as 
zero,  the  DIR  degree  of  freedom  at  all  nodes  not  constrained  in  the  DIR 
direction  will  be  retained  as  masters. 

Master  degree  of  freedom  definitions  are  the  last  input  data  required 
by  the  program.  After  reading  in  this  data  set  the  program  will  print  a 
cross  reference  of  degree  of  freedom  index  numbers  along  with  the  total 
number  of  master,  slave  and  constrained  degrees  of  freedom.  The  mass  and 
stiffness  matrices  are  generated  such  that  their  rows  and  columns  are 
arranged  in  the  order  specified  by  the  cross  reference. 

TRADYR  -  PROGRAM  DESCRIPTION  AND  INPUT  FORMAT 

TRADYR  is  a  special  purpose  program  that  numerically  integrates  the 
equations  of  motion  for  a  plate  and  impactor  during  a  low  velocity  impact 
event.  The  program  assembles  the  mass  and  stiffness  matrices  for  the  plate 
using  the  mass  and  stiffness  matrices  for  cells  or  substructures  generated 
by  LAMCEL  or  LAMSUB.  The  equations  of  motion  for  the  system  are  generated 
using  these  assembled  matrices  and  the  parameters  which  describe  the 
impactor . 

The  complete  set  of  simultaneous  equations  is  reduced  to  a  set  of 
equations  expressed  in  terms  of  a  smaller  number  of  'governor'  degrees  of 
freedom  which  are  defined  by  the  user.  The  reduced  equations  of  motion  are 
then  integrated  numerically  for  the  dynamic  response  of  the  system. 

Input  data  is  read  from  an  ASCII  file  FILE01.DAT  and  from  binary  files 
written  by  LAMCEL  or  LAMSUB  unless  a  restart  analysis  is  being  performed. 

For  restart  analyses  input  is  read  from  FILE01.DAT  and  from  the  binary  files 
FILE03.DAT  and  FILE04.DAT.  These  restart  files  are  generated  at  the  end  of 
the  original  analysis  and  updated  at  the  end  of  each  restart  analysis. 

Output  data  is  written  to  the  ASCII  file  FILE02.DAT,  and  optionally,  to 
a  binary  file  specified  by  the  user.  The  binary  file  contains  the 


A*  21 


displacement  vector  for  each  substructure  at  each  print  atep.  This  file  is 
necessary  for  subsequent  post  processing  in  either  LAMCEL  or  LAMSUB. 

Details  of  the  required  input  data  are  given  in  the  following  sections. 

A.  Control 

TITLE  78  character  alphanumeric  title. 

NCHECK , NRSTRT , NDOF  The  data  check  flag,  analysis 

restart  flag  and  the  number  of 
degrees  of  freedom  per  node. 

NCHECK  -  0  (default),  data  check  only. 

-  1,  analysis 

NRSTRT  -  0  (default),  new  analysis. 

-  1,  restart  old  analysis. 

0  <  NDOF  <  13 

INTALG , NTSTEP , NPRINT , NPSAVE  The  numerical  integration 

algorithm,  the  number  of  time 
steps,  the  number  of  time  steps 
between  solution  printouts,  and 
the  solution  save  flag. 

0  <  INTALC  <  3 

NPSAVE  -  0  (default),  save  solutions. 

-  1 ,  do  not  save  solutions. 

1DPRNT , IVPRNT , IAPRNT  The  displacement,  velocity  and 

acceleration  print  flags. 

IxPRNT  -  0  (default),  print  solution 

-  1,  do  not  print  solution 

DELTA  The  time  increment. 
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NAME1 


The  name  of  the  file  in  which 
the  solutions  will  be  stored. 

This  entry  is  only  included  if 
NPSAVE  -  0. 

For  restart  analyses  (NRSTRT  -  1)  no  other  input  is  required.  However, 
the  analysis  restart  files  FILE03.DAT  and  FILE04.DAT  written  during  the 
previous  analysis  must  be  provided.  At  the  end  of  the  analysis  new  restart 
files  will  be  written  over  those  originally  provided.  If  the  analysis  is 
not  a  restart,  the  program  will  continue  to  look  for  input  on  F1LE01.DAT. 

The  parameter  NDOF  is  not  needed  for  restart  analyses  and  will  be  ignored  if 
included . 

The  program  will  perform  NTSTEP  time  steps  in  the  numerical  integration 
of  the  equations  of  motion  expanding  and  printing  every  NPR1NT  th  step.  A 
time  step  which  will  be  expanded  and  printed  is  referred  to  as  a  print  step. 
If  NPSAVE  has  been  set  to  zero,  the  displacement  solution  for  every  print 
step  will  be  saved  in  binary  on  the  file  specified  by  the  parameter  NAME1 . 

If  NPSAVE  has  been  set  to  1,  these  solutions  will  not  be  saved  and  the 
parameter  NAME1  is  not  Input. 

The  amount  of  information  printed  at  each  print  step  depends  upon  the 
values  of  the  print  flags  IDPRNT,  IVPRNT,  and  IAPRNT.  These  flags 
correspond  to  the  printing  of  the  expanded  displacement,  velocity,  and 
acceleration  solutions  respectively.  When  set  to  zero,  the  corresponding 
data  set  will  be  included  in  the  printout.  When  set  to  1,  the  data  set  will 
be  omitted. 

B.  Model  Definition 

NE  The  number  of  substructures  in 

the  model.  The  program  will 
look  for  NE  sets  of  substructure 
input  after  this  entry. 

NAME2  The  name  of  the  file  from  which 

the  mass  and  stiffness  matrices 
will  be  read. 
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NNE 


The  number  of  nodes  in  this 
substructure.  The  program  will 
look  for  NNE  node  numbers 
following  this  entry,  arranged 
in  groups  of  ten  per  line. 

(NODES (I) ,1-1, NNE)  The  nodes  defining  this 

substructure . 

The  dynamic  model  consists  of  a  number  of  substructures  or  cells  for 
which  the  mass  and  stiffness  have  been  previously  defined  and  stored.  This 
section  of  TRADYR  assembles  the  mass  and  stiffness  matrices  for  these 
substructures  into  a  model  of  the  complete  structure.  The  degrees  of 
freedom  in  this  model  are  the  master  degrees  of  freedom  for  the  individual 
substructures . 

To  ease  the  effort  in  defining  the  substructures  (cells)  that  comprise 
a  model,  each  substructure  is  assumed  to  be  defined  by  a  given  number  of 
nodes  each  having  NDOF  degrees  of  freedom.  These  degrees  of  freedom  are  the 
master  degrees  of  freedom,  specified  for  the  substructure  in  either  LAMCEL 
or  LAMSUB.  This  method  assumes  that  the  same  degrees  of  freedom  were 
retained  for  all  nodes  at  which  masters  were  defined.  If  this  is  not  the 
case  NDOF  should  be  set  to  1  and  each  'node'  on  the  substructure  is  actually 
a  degree  of  freedom. 

LAMCEL  and  LAMSUB  assign  index  numbers  to  master  degrees  of  freedom 
based  on  the  order  in  which  they  were  input,  and  arrange  the  rows  and 
columns  of  the  reduced  mass  and  stiffness  matrices  to  reflect  this  order. 

The  definition  of  the  nodes  which  define  the  substructure  in  TRADYR  must 
also  use  this  same  order. 

C.  Constraints 

NZDF  The  number  of  constrained 

(fixed)  displacements.  The 
program  will  look  for  NZDF 
constraint  definitions  following 
this  entry. 
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NODE, DIR 


NCS 


NIS.DIR 


NODE(I) ,1-l.NIS 


NCE 


NTE 


NODE, DIR, VALUE 


The  node  which  is  constrained, 
and  the  direction  in  which  it  is 
constrained. 

The  number  of  coupled  sets,  the 
program  will  look  for  NCS  sets 
of  coupled  set  input  following 
this  entry. 

The  number  of  nodes  in  the  set, 
and  the  direction  in  which  they 
are  coupled. 

Index  numbers  of  the  coupled 
nodes ,  input  in  groups  of  10  per 
line . 

The  number  of  constraint 
equations.  The  program  will 
look  for  NCE  sets  of  constraint 
equation  input  following  this 
entry. 

The  number  of  terms  in  the 
equation.  The  program  will  look 
for  NTE  constraint  equation 
terms  following  this  entry. 

Constraint  equation  term. 


Constraint  definitions  in  TRADYR  are  interpreted  in  the  same  manner  as 
constraint  definitions  in  LAMCEL  or  LAMSUB.  The  user  is  referred  to  Section 
E  of  the  LMACEL  and  LAMSUB  program  descriptions  for  details  of  these 
definitions.  The  only  exception  being  the  set  of  valid  direction 
specification  labels. 
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Direction  specifications  (DIR)  are  three  character  labels  of  the  form 
Una,  where  n  and  m  are  integers  whose  values  depend  upon  the  value  of  NDOF. 
If  NDOF  is  six  or  twelve  the  direction  convention  defined  in  Section  E  for 
two-  or  three-dimensional  analyses  is  maintained.  If  NDOF  is  any  other 
value,  m  is  a  two  character  integer  ranging  from  01  to  12.  Only  the  first 
NDOF  integers  are  valid  labels,  and  a  zero  must  preceed  integers  less  than 
ten.  Directions  are  always  specified  in  nodal  coordinates  which  may  not  be 
parallel  to  the  global  coordinates. 

D.  Governor  Degrees  of  Freedom 

NGDF  The  number  of  governor  degrees  of 

freedom.  The  program  will  look 
for  NGDF  governor  degree  of 
freedom  specifications  following 
this  entry. 

NODE, DIR  The  node  number  and  direction 

that  define  the  governor  degree 
of  freedom. 

The  direction  parameter  DIR  may  be  any  valid  label  as  defined  in  the 
constraint  section  (Section  C) .  In  addition,  DIR  may  be  assigned  the  label 
ATT,  in  which  case  all  unconstrained  degrees  of  freedom  at  the  specified 
node  will  we  retained  as  governor  degrees  of  freedom.  If  NODE  is  input  as 
zero,  the  DIR  degree  of  freedom  at  all  nodes  not  constrained  in  the  DIR 
direction  will  be  retained  as  governors. 

E.  Impact  Parameters 

NODE, DIR  The  node  at  which  the  impact 

occurs  and  the  direction  of  the 
impact . 

MIMP,K,N,M, I ,C  The  mass  of  the  impactor  and  the 

five  spring  constants. 
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The  Initial  position  and 
velocity  of  the  impactor. 
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